Non-Equilibrium Electron Transport in Two-Dimensional Nano-Structures Modeled 
by Green's Functions and the Finite-Element Method 



O 

o 
< 



I 

O 

o 



> 
o 
oo 

o 

cn 

o 



I 

O 

o 



p. Havu\ V. Havu^, M. J. Puska\ and R. M. Nieminen^ 
1) Laboratory of Physics, Helsinki University of Technology, P.O. Box 1100, FIN-02015 HUT, Finland 
2) Institute of Mathematics, Helsinki University of Technology, P.O. Box 1100, FIN-02015 HUT, Finland 

We use the effective-mass approximation and the density-functional theory with the local-density 
approximation for modeling two-dimensional nano-structures connected phase-coherently to two in- 
finite leads. Using the non-equilibrium Green's function method the electron density and the current 
are calculated under a bias voltage. The problem of solving for the Green's functions numerically is 
formulated using the finite-element method (FEM). The Green's functions have non-refiecting open 
boundary conditions to take care of the infinite size of the system. We show how these boundary 
conditions are formulated in the FEM. The scheme is tested by calculating transmission probabilities 
for simple model potentials. The potential of the scheme is demonstrated by determining non-linear 
current-voltage behaviors of resonant tunneling structures. 
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I. INTRODUCTION 



Two-dimensional (2D) nanodevices are structures in 
which electrons move in a restricted nanometer-size area. 
The phase-coherence length of electrons is of the order of 
the dimensions of the device. Electron transport through 
nanodevices cannot be modeled using the traditional de- 
scription based on diffusion or Boltzmann equations. One 
has to use a method which takes the quantum-mechanical 
character of the carriers, e.g. quantum interference, ex- 
plicitly into account [1]. 

Nanodevices are fabricated using semiconductor- 
heterostructure techniques. A layer of semiconductor 
(e.g. AlGaAs) is grown on top of another semiconductor 
(GaAs) with molecular-beam epitaxy. The two semicon- 
ductors have different band gaps so that electrons accu- 
mulate in the potential well at the semiconductor inter- 
face and form a 2D electron gas. Above of the semicon- 
ductor layer metallic gates are fabricated. Applying volt- 
age on them the electron motion can also be restricted in 
the horizontal direction and nanodevices, such as quan- 
tum point contacts and quantum dots, are created. 

The quantum-mechanical modeling of 2D nanostuc- 
tures is usually based on the effective-mass approxima- 
tion. For the ground-state carrier distribution one can 
employ, for example, Monte Carlo-methods [2] or density- 
functional theory (DFT) [3] . The description of isolated 
structures is rather straightforward because the system 
is finite and all the electron states can be calculated. 
Often the nanodevice is connected to a measuring sys- 
tem by leads and the current through the system is mea- 
sured. If the connection is weak the nanostucture can 
still be approximated as an isolated system, but in the 
case of strong coupling the combined nanostructure-leads 
system has to be described. In this case the leads can 
have a considerable effect on the electronic structure of 
the nanodevice. The electronic structure of this kind of 
open system can be obtained using DFT by calculating 



the wave functions in the scattering formalism using the 
Lippmann-Schwinger equation [4]. The method also re- 
lates to the conductance of the system in the limit of zero 
bias. Another possibility is to use DFT in combined with 
the non-equilibrium Green's function (NEGF) method 
[5] . In this scheme the wave functions are not calculated 
explicitly in the device region. The NEGF-approach also 
enables the addition of a bias voltage between the leads 
and the calculation of the current through the system 
also in the non-equilibrium state. 

The electronic-structure calculations using the Green's 
functions demand extensive computer resources. There- 
fore the numerical method for the Green's function im- 
plementation has to be chosen carefully. There is a wide 
range of different numerical methods available today for 
electronic structure calculations, e.g. the finite-difference 
method [6], the linear combinations of atomic orbitals 
(LCAO) method, the wavelet method [8], and the plane- 
wave method [9] among the most popular ones. Previ- 
ously, the Green's function method coupled to DFT has 
been used in nanostucture calculations employing atomic 
orbitals [10,11], localized optimized orbitals in real space 
[12], Gaussian orbitals [13] or wavelets [14] as basis func- 
tions. 

In the present work we have adopted the finite-element 
method (FEM) to study 2D nanostructures within the 
effective-mass theory and using the DFT-NEGF scheme. 
Previously, electronic structure calculations the FEM has 
been used in, for example, in Refs. [15-19]. The main ad- 
vantages gained by the FEM in the present context are 
the possibility to control the accuracy of the approxima- 
tion via mesh refinements, the ability to simulate easily 
different geometrical configurations of the system and the 
ease in the treatment of the boundary conditions. More- 
over, the evaluation of the basis functions is fast and 
the ensuing sparse linear systems allow the use of fast 
sparse solvers. In practice, we have chosen to use piece- 
wise polynomials as basis functions. The polynomials are 
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very fast and stable to evaluate in any computational en- 
vironment. The approximation properties of the polyno- 
mials are well-known and several error bounds are avail- 
able [20] . In the FEM the open boundary conditions are 
easier to implement than in the finite-difference method 
[1] and in the basis set methods [10,11,14] in which they 
are derived by first writing down the infinite discretiza- 
tion matrix and then cutting out the central area from 
it. In the FEM these boundary conditions are written in 
a simpler and more intuitive way as will be shown in this 
work. 

We use use effective atomic units which are derived by 
putting the fundamental constants e — h = rUe = 1, and 
the material constants, the effective electron mass and 
the dielectric constant m* = e = 1 respectively. The ef- 
fective atomic units are transformed to the usual atomic 
units using the relations 

Length: 1 aS = 1-^ao w ^ 0.529177 x lO^i" m 
Energy: 1 Ha* = l^^Ha 27.2116 eV 

Current: la.u.* = l|^a.u. « 2^6.6231 mA. 

The organization of the present paper is as follows. 
In Sec. II. we present our 2D nanostucture model and 
explain how the Green's functions are used in the elec- 
tronic structure and current calculations. In Sec. Ill we 
formulate the solution of the Green's functions within 
the FEM. Finally, in Sec. IV we deal with our test cases, 
which include confining well and bottle-neck model po- 
tentials and double-wall barrier systems. Sec. V contains 
the conclusions. 



II. MODEL AND GREEN'S FUNCTION 
FORMULATION 

A. The model for two-dimensional nanostructures 

In real nanodevices electrons of the 2D electron gas 
are in a potential well at the interface between two semi- 
conductors. The electron density in the well is neutral- 
ized by a positively-ionized donor layer separated from 
the potential well. The lateral confinement of electrons 
is obtained by gate voltages. Electrons are in practice 
in the ground-state with respect to the motion perpen- 
dicular to the interface. Therefore our model is strictly 
two-dimensional. 




FIG. 1. Model nanostructure between two infinite leads. 

A schematic sketch of the model is in Fig. 1. It shows 
the region of interest between two semi- infinite leads. 
The potential profile is a combination of interactions be- 
tween electrons and the positive constant background 
charge (jellium), and the external potential caused by the 
gate voltages. Thus, the layer of ionized donors and the 
2D electron layer coincide in our model. In many models 
the potential profile is approximated using a harmonic 
potential profile [3,21]. In our model this approximation 
cannot be used, because we solve for the electrostatic 
potential of an infinite system requiring that the system 
is charge neutral. In order to keep the model simple the 
confinement of the electrons is established by shaping the 
background charge and, optionally, by external potentials 
in certain regions of the system. 

We divide the infinite system to three separate areas 
as shown in Fig. 1, the central area fl, the left region 
57i and the right region fi^j. We denote the boundary 
between the regions and flu as d^R and between the 
regions and as dflR. The Green's functions are cal- 
culated in the region fl. dflL and OQr, are non-reflecting 
open boundaries. On the other two boundaries, dflpi^p2 
which are far enough from the important device region 
the potential is assumed to be infinite, so that the Green's 
functions vanish there. 

We solve for the self-consistent electron structure of 
the system iteratively. The electron density is calcu- 
lated from the Green's functions. The effective potential 
is calculated from the electron density as usual in the 
DFT within the local-density approximation (LDA). Af- 
ter mixing the new effective potential with potential from 
the previous iteration the electron density is recalculated. 
The loop is repeated until convergence is achieved. 

The effective potential has four terms 

Vef f^Vc + V^c + Vuas + Vgate , (1) 

where Vc and Vxc are the Coulomb and the exchange- 
correlation potentials arising from the charge distribu- 
tions, respectively. The calculation of Vc is discussed be- 
low in more detail. For Vxc, we use the recent 2D-LDA 
functional by Attacalite et al. [22,23]. 

Vbias takes care of the boundary conditions under the 
bias voltage [5] . The total electrostatic potential has dif- 
ferent levels in the right and left leads. This introduces 
Vbias as a linear ramp potential over fi. In the regions ^II 
and Qr, Vcff is calculated as a potential of the infinite 
(jellium) wire. Then Ve// is also continuous if f2 is large 
enough. The ensuing energy scheme is shown in Fig. 2. 
Also the Fermi levels in the right and left leads differ by 
the applied bias voltage AVbias- Vgate is an external gate 
potential. Using gate voltages it is possible to increase or 
decrease the potential in certain regions, for example to 
increase the potential walls and to decrease the potential 
wells of a bare jellium system. 
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FIG. 2. Effective potentials and Fermi levels under the bias 
voltage. 

Below wc use a notation in which a point inside the 
two-dimensional region fl is denoted by r and a point 
outside the region O in region flu or Ql hy re- A point 
on the boundary dQ,L is Tl and a point on OSIr is vr. 



B. Green's functions in electronic-structure 
calculations 



We use Green's functions in calculating the electronic 
structure and the current under an external bias voltage. 
The theory is explained in more detail in Refs. [1] and 
[5]. The electron density is calculated from the Green's 
function . In order to obtain one has to solve first 
for the retarded Green's function C from 

{L0-H{r))G''{r,r';w)=5{r-r'), (2) 

where cu is the electron energy and H is the DFT Hamil- 
tonian the system, 

H{r) = -lv' + Veffir). (3) 

In this case r is a two-dimensional variable. Its compo- 
nents along and perpendicular to the leads are x and y, 
respectively. C is zero on the boundaries parallel to the 
leads (sec Fig. 1). If cj is smaller than the bottom of the 
potential 14// in the lead Eq. (2) gives exponentially de- 
caying solutions there. Otherwise the solution oscillates 
with a constant amplitude to the infinity. The form of 
G'~{r, r') in a uniform jellium wire is shown in Fig 3. The 
real part has a pole at r = r', while the imaginary part 
behaves smoothly everywhere. This is why the imaginary 
part is much easier to approximate numerically than the 
real part. 



a) 
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FIG. 3. Real (a) and imaginary (b) parts of the Green's 
function G^{r,r') for a uniform jellium wire, r = {x,y) and 
r' = (21.6, 15.4) (the position of the pole). 

In equilibrium, when the Fermi functions in and 
flu are identical, fhiw) = /fl(w), we obtain 

G<(r,r» = 2fL,R{w)G-{r,r';w)- (4) 

This equation is also valid under a bias voltage at ener- 
gies uj for which /l(w) = fniuj) (in practice, fL/R = 1 
for those energies). If Eq.(4) is not applicate, has to 
be calculated in a more complicated way. Eq. (2) can be 
reformulated using the so-called retarded self-energies of 
the leads, and S£, as 

(a; - i?o - Si (a;) - S5^(a;))G'-(r, r'; w) = 5{r - r'). 

(5) 

Above, Hi) is the Hamilton operator for the isolated cen- 
tral area f2. In practice, T^i^/r can be calculated from 
the boundary conditions for the Green's functions at 
d^L/R. ^L/R are functions with non-zero values only 
at the boundaries dfl^/R. Next we define the functions 
^L/R as 

iVL = - Si = 2z3(S£), 

zTr = S^, - S?, = 2*3(S^j). ^ ' 

'^L/R. self-energies for the advanced Green's func- 

tion C = {G^y* . One can then write the electron density 
as the sum of the electron flows from the leads to the re- 
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gion CI, using 
G<(r,r';a;) = 

-ifniu) [ [ G^r,rR;u>)TR{rR,r'j,;u) 

X G''{r'j^,r'-uj)drRdr'j^ 

r r (7) 

JdQL JdQL 

xG%r'L,r';co)drLdr'L, 

where Jr/l are the Fermi functions in the right and left 
leads. This equation has to be used in nonequilibrium 
situations when Jr^ /l- 

Eq. (7) corresponds to the electron density due to the 
states extending to infinity in the leads. Eq. (4) includes 
also the electron density of possible bound states, which 
are localized near O and decay exponentially in the leads. 
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FIG. 4. Integration path used in Eq. (8). 

In order to calculate total electron density we integrate 
over the electron energy u 

P(r) = ^ j_ S(G<(r,r;a;))cL;. (8) 

We use both equations (4) and (7) in this integration. 
Eq. (4) is analytic in the upper half of the imaginary uj- 
plane whereas Eq. (7) has poles below and above the real 
ui axis. Thus, using Eq. (4) it is possible the transfer the 
integral path from the real axis to the complex plane. 
Our integration path is shown in Fig 4. The first part 
is a semi-circle Cj in the complex w-plane using Eq. (4) 
and it takes care of the possible bound states below the 
energy bands of the leads. The rest of the integration, 
C//, is close to the real axis and there Eq. (7) is used. 
On the semi-circle only few integration points are needed 
because the rapid variations of are smeared out when 
the integration leaves the real axis. This is specially use- 
ful for the bound states, which give rise to sharp peaks 
near the real axis. 

Computationally, it is faster to solve for G^ from 
Eq. (7) than from Eq. (4). Eq. (4) results in the inversion 
of the entire matrix, because one needs G'^{r, r') in all the 



discretion points of Q. Electron density in Eq. (8) is the 
calculated using the diagonal entries of the imaginary 
part, 9(G''(r, r)). Inversion of the matrix using direct 
sparse routines from HSL [7] occurs as follows. First one 
performs the symbolic analysis and factorization to pro- 
duce an ordering that reduces the fill-in. After that a 
numerical factorization with pivoting is performed pro- 
ducing the Cholesky factor of the matrix. The a set of lin- 
ear equations with different right-hand sides are solved. 
The number of equation is equal to the dimension of the 
matrix. Eq. (7) needs only the Green's functions G^{r, r') 
for r' = Ti/jf on the boundaries dflL/R. This means that 
after factorization one has to solve for a set of only as 
many linear equations as there are discretization points 
on dflL/R- 

For 2D systems the use of Eq. (4) is justified because 
the analytic continuation of the integrand reduces the 
number of points needed in the numerical integration of 
Eq. (8) and because the discretization error is smaller for 
Eq. (4) than for Eq. (7). Namely, only the imaginary 
part of C is used in Eq. (8) so that the pole of 5R(G'') 
does not cause any major numerical problems if Eq. (4) 
is used. 



C. Electric Current 

The electric current is also calculated using the Green's 
functions. The electron tunneling probability through 
the central region is obtained from 

T{to)=[ [II TL{rL,r^;u;)G^r'L,rR;u;) 

JdQL JdSlL JdQ.[i JoSIr 

xVRiTR, t'r, uj)G\r'R, rL\uj) dvL dr'^ dvR dr'R, (9) 

and the total current is calculated integrating over the 
energy, u and taking care of the electron occupations in 
both leads. In the effective atomic units the result is 

I= - r T{ij) {fL{oj) - fR{w)) dw. (10) 

7^ ^-oo 



III. FINITE-ELEMENT METHOD FOR SOLVING 
GREEN'S FUNCTIONS 

A. Variational formulation 

The most demanding computational task is to find the 
Green's function at different energies as presented above. 
To this end, we first divide the domain of the problem 
into two disjoint parts, the computational domain fl and 
the exterior domain 17^. Only the computational domain 
is discretized whereas the exterior is taken care of by 
the corresponding Green's function (see below Gh. HIE). 
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First, wc cast Eq. (2) into a variational, or weak, for- 
mulation for the domain f2. During the derivation we 
frequently make use of the Green's formula 



/ Vu • S/v dr = I ds — I uS/^v dr, 

Jn Jan dn Jq 



(11) 



for two arbitrary, sufficiently smoothly behaving func- 
tions u and V. Above, n denotes the outward normal 
of f2, and the line integration is taken in the counter- 
clockwise direction around the 2D area fl. 

To proceed, we multiply Eq. (2) by a sufficiently 
smooth function v and integrate the resulting identity 
over f2 giving 

/ v{r) [lj - H{r)] Cir, r'; u) dr 
Jn 

+ [a;-ye//(r)]G'-(r,r»}dr ^^^^ 

= / v{r)5{r — r') dr 
Jn 

= v{r'), 

The use of the Green's formula of Eq. (11) gives 



/ v{r))-V'^G''{r,r';uj)dr 
Jn 2 

= - / Vv{r) ■\vG'{r,r']uj)dr 
Jn 2 

15G'-(rL,r';w) 



+ 



drL 



/ v{rL)-- „ 
ion,, 2 duL 

JanR 2 dnR 



(13) 



Thus, the original problem of Eq. (2) is equivalent to the 
formulation 



1, 



- Vv(r) • -VG"'(r,r';a;) 
n >■ 2 



+ v{r) [lj - Vef f (r)] G"^ (r, r» } dr 

f ldG''{rL,r';uj) /..n 

+ / 7^ -v{rL)drL y^^) 

ion. 2 driL 



an. 



ldG^irR,r';uj] 
2 driR 



v{rR) drt 



-.v{r') 



for any sufficiently smooth function v. 

In order to obtain a solvable system, the boundary 
conditions must be supplied at the boundaries OQl and 
dflR. For conciseness we discuss only the case of O^Il, 
the other case OVIr being similar. Gonsider the exterior 
problem 



{u - H{re))ge{re,r'^;uj) = Sire - r'J, r'^ € CIl 

9e{re,K;uj) =0, re&dflL, (15) 

for the Green's functions Qe of the semi-infinite lead. It 
follows that any sufficiently smooth function u can be 
written in the form 



''(K) = / u{re)S{fe -K)dre 

J nr. 



u{re)[oJ - H{re)]ge{re,r'^;oj)dre 

^*(^e){^VV(r-e,<;a;) (1^) 

+ ['^-yeff{r)]geire,K;to)^dre 

for Tg G flL- Using the Green's formula (11) for the ex- 
terior domain Q,l twice we can write 

/ u{re)l:V'^ge{re,ri;ui)dre 
JnL ^ 

= - I ^Vu{re) ■'Vge{re,K;i^)dre 
JnL ^ 



an, 



^u[r^) dr^ 



- [ ^ffe(r-e,r-e;a;)V^M(re)(ire 

J nr. ^ 



(17) 



on. 2"^'^^) a< 



an^ 2 dn^ 



T^9e{r'L,r'e,i^)dr'i^, 



so that 



"(^e) = / ge{re,r'^;uj)[uj - H{re)\u{re)dre 
JnL 

1 ^99e{r'L,r'^;uj) 



+ 



JanL ^ 



-i 



(^^L (18) 



^^''^''-^9e{r'„rUco)dr', 



anL 2 dn'L 



We assume that u = G'' is the solution to the homoge- 
neous problem (oj — H{re))G'^{re,r'-,oj) = for G fiz,. 
Since (?e = on d^L we have by Eq. (18) 

G'-(r;,r» = / iG'-(rLr';a;)^^^%^<, 
JanL ^ '^^L 



r' e nL. 



(19) 



Now the representation formula (19) can be used to sup- 
ply the boundary condition to Eq. (14). Differentiating 
Eq. (19) with respect to and letting — > ri, e dQ,L 
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we obtain the term corresponding to the left boundary 
dflL in Eq. (14) as 



L 



ldG^{rL,r';uj) 
2 driL 



v{rL) drL 



dviL Jen 



driLdn'j 



xv{rL) dr'^drL 



(20) 



Here we have derived the variational form for the self- 
energy operator Ej,. It includes line integrals over the 
boundary OQl together with a trace mapping from func- 
tions on n to the functions on d^L- The function S£ in 
Eq. (6) is given by 



1 d'^ge{r'j^,rL;(^) 

2 driLdn'r 



(21) 



with zero extension outside the boundary 917^. 

The mapping generated above by Eq. (20) is called the 
Dirichlet-to-Neumann mapping since in general it maps 

the Dirichlct datum u of a solution to a partial differen- 
tial equation to the corresponding Neumann datum |^ . 



B. Finite-element discretization 

To obtain a numerical approximation for the Green's 
function G"" in the computational domain O wc select a 
finite-dimensional space Sh defined on f7 and project our 
problem of Eq. (14) into Sh by solving for GJ^ e Sh such 
that 



IS 



1 



2 

+ 



WGl{r,r';u:)-Wvh{r) 

[uj - V^ff{r)\Gl{r,r';uj)vh{r)^ dr (22) 

+ <tLGl,Vh> + <tRGl,Vh> 
= Vh{r') 

for every Vh G Sh [28] . A matrix equation is obtained by 
selecting a basis {4>i}fLi for Sh and expanding G^ in the 
basis, 

GUry)=J2alM^)4>jir'). (23) 
Selecting Vh = 4>k in Eq. (22) we obtain 

N 

+ [oj- Veff{r)]cl)i{r)Mr)] dr (24) 

-I- < t,L(l)i, (j)k> + < t,R(j)i, 0fe > I 



Denoting 

aik = J^(-l^Mr)-^Mr) 

+ [oj-Veff{r)]<Pi{r)Mr))dr ^^5) 

-I- < Y.L(l)i,(l>k > + < ^R(l>i,(l>k >, 

and exploiting the symmetry of the coefficients gij we see 
that gij 's are the entries in the inverse of the matrix the 
given by Eq. (25). 

We connect S/,/^ to the discretized forms as 



^L/B.,t,j — < ^L/R4'ij4>j > ■ 



(26) 



Further, let us denote 

Gl = Y.f^kiMr)Mr'), (27) 



k.l 



and 



with 



t^/R = 29(S2/^), 



L/R,ij 



L/R,ji 



(28) 



(29) 



since r^/7j is symmetric. Now, for example, the electron 
tunneling probability of Eq. (9) can be written in the 
discretized form as 



n-)= I I I j 

i,j,k,i=i''^^^^ "'ant Jdnn Jdnn 



TL{rL,r'L)9liHr'L)'p3{rR) 
X VR{rR, t'r) gli (pkivR) M^l) 
X drLdrLdvRdr'R, 

N ^ ^ (30) 

= X] < ^^^i' > Slj < ^R(t>k, (t>3 > 9kl 
i,j,k,l=l 
N 

= XI ^^''^ 3ij ^R,jk gti ■ 

i,j,k,l=l 

C. Finite-element basis 

So far we have not touched the subject of selecting the 
basis functions (pi in Sec. B above and thus the space Sh- 
in principle, we could select any computable set 
but adhere to a traditional choice in the finite-element 
practice, namely to the set of piecewise polynomial func- 
tions. The basis functions arc constructed as follows. As- 
sume that il is partitioned into a simple mesh of N nodes 
and M polygons Tj conforming to the usual requirements 
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imposed on a finite-element mesh. These polygons can 
have a variety of shapes but the simplest choice of tri- 
angles in two (and tetrahedral in three) dimensions will 
serve our purposes. We choose the basis functions (j^i to 
be element- wise linear functions that have the value one 
in a single node of the mesh and zero in other nodes (see 
Fig. 5). The corresponding finite-element space Sh is 

N 

Sh={vh = ^Ci(j)i\ci eC} ^^^^ 

={vhGC{n)\vhiT,€Vi{Ti)}, 

where C (f2) denotes the set of continuous functions in f2 
and ■Pi(Ti) is the set of polynomials of degree one in the 
polygon Ti. 

An element-wise polynomial basis has several advan- 
tages. First, polynomials are fast to evaluate and they 
can be integrated exactly on a suitable reference element. 
Second, the piecewise nature of the functions ensures that 
the matrix (aij)ij=i is very sparse. Third, the accuracy 
of the discretization can be controlled via mesh refine- 
ments and coarsening. 

The piecewise nature of the basis functions gives rise 
to a sparse matrix. Due to recent developments in linear 
algebra there are fast direct solvers [24] (also parallel) 
[25,26] for sparse systems arising from discretization of 
partial differential equations. Since we must solve for all 
the coefficients gij of the approximate Green's function 
gh we are faced with the problem of solving N linear 
systems with different right-hand sides. This kind of set- 
ting is favorable to direct methods over iterative ones. 
Nevertheless, the computation itself is a time-consuming 
procedure and cannot be substantially accelerated with 
the techniques known today. 




FIG. 5. A linear basis function (p. The function is one in a 
given mesh node and descends linearly to zero in the adjacent 
nodes. 



D. Mesh generation 

An important property affecting the quality of the 
finite-element approximation is the underlying mesh and 
especially the shape and the size of individual elements. 
Several techniques for mesh generation in two and three 
dimensions are available. All the techniques have in com- 
mon that they try to produce meshes with elements of 



desired local size and high quality. There are also sev- 
eral indicators for evaluating the quality of the shape of 
a single element. Perhaps the most common is to reqiiire 
that there are no large angles in the element. Typically, 
the larger the maximal angle of an element is, the worse 
the resulting approximation will be. 

In this work we use Delaunay meshes [27] for trian- 
gular elements in two-dimensional problems. They are 
known to be very robust in producing high-quality trian- 
gular meshes for different shapes of domains. A Delaunay 
mesh can be characterized as follows. A mesh consisting 
of N nodes and M triangular (or tetrahedral) elements 
satisfies the Delaunay criterion if the circumscribe Cj of 
a triangle (or tetrahedron) Tj of the mesh contains no 
nodes of the mesh. Meshes satisfying the Delaunay cri- 
terion are called Delaunay meshes. 

It can be shown that for a given set of points in a plane 
a Delaunay triangulation always exists and is even unique 
with a minor assumption on the placement of the nodes. 
Furthermore, among all triangulations of the nodes, the 
Delaunay triangulation maximizes the minimum angle 
present in the triangulation. The max-min property can 
be usually considered as a guarantee of high-quality ele- 
ments. 

Unfortunately the Delaunay criterion is not sufficient 
for a high quality tetrahedral mesh in three dimensions. 
This is due to the presence of "slivers" in Delamiay 
meshes. These elements can have very large angles deteri- 
orating the approximation capabilities, and yet they sat- 
isfy the Delaunay property. Therefore alternative tech- 
niques must be sought for when producing meshes in 
three dimensions. Typical approaches use a mixture of 
different methods, e.g. octree methods, advancing front 
methods, and Delaunay methods. 

However, it should be noted that the quality of the 
resulting mesh produced by a mesh generation algo- 
rithm depends heavily on the shape of the domain to be 
meshed. Very simple domains such as cubes and other 
rectangular domains are usually well treated by virtually 
any method, whereas more complicated domains having 
holes and cuts need more attention. 



E. Exterior Green's function 



The exterior Green's function for the semi-infinite 
leads can be calculated numerically as the surface Green's 
function of a periodic system [14]. In the present work 
the potential is uniform in the leads along the lead axis. 
Therefore we can solve for the isolated Green's function 
using the analytic one-dimensional solution along the 
lead and the numerical transverse wave functions Xm{y) 
[1]. The ensuing exterior Green's function for the quasi- 
two-dimensional semi-infinite wire is 
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1=1 ^ '(32) 

where Xm(2/)'s are solutions to the Kohn-Sham equation 

. i - Ve// (y) j Xm (2/) = emXm {v) , (33) 



with 



km = 1/2(0; - €n 



(34) 



We solve Eq. (33) using self-consistency iterations for the 
electron density and the potential profile Veff{y)- As 
explained before we use a model in which the positive 
charge forms a thin wire and the electron wave functions 
spread out of this charge. The effective potential V^ff 
consists only of Vxc and Vc, and no external potential is 
applied. In practice the summation in Eq. (32) is trun- 
cated typically after a few tens of states so that the results 
are well- converged. 

The charge densities resulting from this calculation are 
used in the boundary conditions when calculating the 
Coulomb potential of the nanosystem. The total charge 
per unit length is zero in an infinite wire, but there are 
local variations in the charge density in the transverse 
direction. As an example, we show in Fig. 6 the effective 
potential and the positive and negative charge densities 
in a case with two transversal modes in the wire. A cut 
perpendicular to the wire axis is shown. 



Coulomb is treated in three dimensions. In this case it 
is not efficient to solve for the three-dimensional Poisson 
equation, but to evaluate the integral 



Above, p is the electron density and pp is the positive 
background charge density. The integral is evaluated by 
integrating basic functions in every element. For ele- 
ments with no pole (r is not inside the clement), the 
integral is evaluated using the Gaussian quadrature rules 
for triangles [29]. Elements which have r in one corner 
are evaluated by making a mapping from the triangle to 
a square in which the pole disappears [30] . 



IV. TEST SYSTEMS 



This section is devoted for testing and demonstrat- 
ing our scheme. First the transmission probability over 
a given potential well and through a given bottle-neck 
potential are determined. The aim of these non-self- 
consistent calculations is to provide, trough the compar- 
ison with the exact results, an idea of the numerical ac- 
curacy of our methods. Thereafter we demonstrate the 
possibilities of the scheme by solving self-consistently the 
electronic structure and the current under a bias voltage 
for different resonant tunneling systems. 
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FIG. 6. Electron density (solid line), positive background 
charge (dotted line) and Vef / (dashed line) for an infinite uni- 
form wire. 



F. Coulomb interactions 

The effective potential is also calculated using the 

FEM and the same mesh as for the Green's functions 
is used. Vxc is simply evaluated in every node point. The 
potential charge densities are two-dimensional but the 



A. Transmission probability over a potential well 



Basic quantum mechanics gives the transmission prob- 
ability over a potential well (see the inset in Fig. 7) as 



1 + 



yo^sin^(V2(u; + Vo)L) 
4uj{uj + Vo) 



(36) 



where Vo and L arc the depth and the length of the well, 
respectively, and u is the electron energy. Our numeri- 
cal approach obeys this result accurately. For example. 
Fig. 7 gives the transmission probability calculated using 
Eqs. (9) and (30) for a narrow wire with a potential well. 
For the energies shown there is only one transverse mode 
in the wire. The good agreement between the numerical 
and analytic results indicates that the FEM mesh is fine 
enough. 
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FIG. 8. Bottle-neck model potential. The potential is con- 
stant inside the leads and in the bottleneck between the leads. 

At the boundaries the potential rises to infinity. The dimen- 
sions are L — H — IOuq and W — 30 Uq. The length of the 
calculation area 5* = 30 . The FEM mesh shown has smaller 
elements near the boundaries diij^/n. 



FIG. 7. Transmission probability over a potential well. The 
solid line corresponds to the analytic solution of Eq. (36) 
and the circles are calculated using the FEM code. In this 
calculation L = 10 Oq, Vb = IHa*, the width of the wire 
W = 2>aQ, and the average distance between the FEM mesh 
nodes h = 0.3 a^. 



B. Transmission probability tiirough a bottle-neck 
potential 



Next we study how the FEM node density affects the 
results. We calculate the electron transmission probabil- 
ity as a function of energy using different FEM meshes. 
Our scattering potential is a bottleneck shown in Fig. 8. 
The electron transmission probability is shown in Fig. 9 as 
a function of the energy. In stepwise jumps in the trans- 
mission probability mean that new transverse modes 
cmerg with increasing energy to. The narrow peaks near 
the beginning of each step correspond to the constructive 
interference of the incident wave with the wave reflected 
twice at the lead-bottleneck boundaries [31]. Increasing 
the energy means making the electron wavelength shorter 
so that more points are needed to describe the wave func- 
tions. Thus, with a fixed element size h it is possible 
to characterize transversal modes up to a certain energy 
only. Thereafter the transmission probability collapses 
due to the a loss of numerical stability. 



L 




In Fig. 9a the size of the elements in each calculation 
is the same throughout the whole calculation area. Ac- 
cording to the two uppermost curves corresponding to 
the FEM node distances h = lag and /i = 2 Aq , we need 
about 4 nodes between the adjacent zero- value lines of 
the electron wave function. This means that the FEM 
node distance of /i = 3 Cq should give a reasonable result 
for the first transversal mode. In contrast, the results 
show large oscillations of the transmission due to dis- 
cretization errors. The reason for this is that the pole of 
the real part of the Green's function is not approximated 
accurately enough. When determining the transmission 
the arguments of the Green's function are on the opposite 
boundaries (Eq. (9)). These Green's function values are 
calculated by solving a linear equation problem in which 
one of the arguments of G"'(r, r') is fixed e.g. on the left 
boundary, d^L and the other argument runs over the 
central region to the right boundary O^Ir. If the FEM 
mesh is not dense enough near the left boundary where 
the pole is a large numerical error propagates to the el- 
ements needed in Eq. (9) [32]. In Fig. 9b the number of 
points at the boundaries dCl l /r is larger than inside the 
calculation area f2. The figure shows that the effects of 
the discretization errors are now strongly reduced at low 
energies, but the transmission probability at high ener- 
gies collapses as fast as in Fig. 9. In conclusion, when one 
wants to describe the transmission probability only up to 
a certain energy value, the optimum way to choose the 
sizes of the elements is to use smaller elements near the 
boundaries dQ^/R than inside the area fl. In this simple 
test system the bottleneck potential is relative wide, but 
if the bottleneck is narrow in comparison with to the rest 
of the wire, it is reasonable to refine the mesh also in the 
neck region. Finally, the above refinement is also needed 
when calculating the electron density in nonequlibrium 
using Eq. (7). The real part of G"'(r, r') is needed be- 
tween a point on the boundary, dflL,R&-TLd an arbitrary 
point in the central region CI. 
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FIG. 9. Electron transmission probability as a function of 
the energy for different FEM meshes, a) All the elements 
in each calculation are of the same size. The FEM node 
distance h = lag (solid line), h = 2aQ (dashed line) and 
h = Sag (dotted line), b) The elements are smaller near 
the boundaries dQ^/ji (see Fig. 8). The minimum distance 
lal and the maximum distance hmax = 2a5 (solid 
line) and hmax = 3a5 (dashed line). 



C. Resonant tunneling through double-barrier 
potential systems 

1. Symmetric harrier system 

In this subsection we demonstrate the potential of our 
scheme by showing results of self-consistent electronic- 
structure calculations for 2D nanostructures under a fi- 
nite bias voltage. We restrict ourselves to zero temper- 
ature calculations. The test system is a double-barrier 
potential structure, a schematic sketch of which is shown 
Fig. 10a. A jellium wire is cut by two vacuum regions and 
additional potential barriers are introduced within them 
in order to adjust the potential and the transmission. We 
consider two special cases. Case A has thinner potential 
walls l'^'" = 1 ^5 than case B for which l'^'" — 1.25 Oq. 



This difference means that the connection to the leads 
differs remarkably its the strength. We make contact 
with real semiconductor systems by converting our re- 
sults from the effective atomic units to the STunits using 
the effective mass of electrons m* = 0.067 and the dielec- 
tric constant e — 12.4 for GaAs. Then Oq = 9.779 nm and 
1 Ha* — 11.8672 meV. The positive background charge 



density 0.2 (ag)" 



2 • 10 nr corresponds to a rea- 



sonable electron density at the GaAs/AlGaAs interface. 
The groundstate electron density of the double-barrier 
system is shown in Fig. 10b, exhibiting Friedel oscilla- 
tions in both leads. The wires are so thin that only one 
transverse mode is occupied. 



a) 



Vw 



Vw 
I 



W 




FIG. 10. Double-barrier potential system, a) The model. 
The gray areas correspond to the positive background charge. 
At the gaps there is an additional potential Vw — 2 Ha* . The 
size of calculation area is 29 x 5(a5)^, the width of the 
background charge W = Sag and length of the quantum dot 
L — 9 oq. Case A has = 1 ao and case B = 1.25 ag. 
The number of FEM nodes used in the calculations is 2105. 
b) The total electron density at zero bias voltage for case A. 



The effective potential along the symmetry axis of the 
double-barrier system at zero bias voltage is shown in 
Fig. 11a. The potential barriers are so small that the 
quantum dot is strongly connected to the leads. When 
we add the bias voltage to the system, the potential of 
right lead increases and that of the left lead decreases. 
The change of Veff for case B is shown in Fig. lib. The 
maximum bias voltage applied is small in comparison to 
the barrier heights. The potential drop occurs between 
the potential walls, not in the leads. This is expected 
because the leads are ballistic, with no scatterers at all. 
At small AVbias values the potential in the quantum dot 
stays at the level of the potential in the left lead. This 
is seen in the upper panel of Fig. lib. When AVbias is 
large enough the potential in the dot rises close to the 
mean value in the leads (see the lower panel). A nearly 
inversion-symmetric potential develops. In case A the 
potential in the quantum dot develops differently. It fol- 
lows mainly the potential level of the right lead for all 
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bias voltages studied. 
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quantum dot stays at the level of the left lead. However, 
when AVbias is large enough the resonance peak enters 
the bias window, the charge redistribution occurs quite 
symmetrically at both barriers and the potential level in 
the quantum dot is in the middle between the left and 
right lead levels. The resonance peak of case A is wider 
than that of case B because the connection to the leads is 
stronger. The wide resonance enters the bias window at 
a low bias value and its position follows the Fermi level 
of the right lead. Then the bias-induced charge redistri- 
bution takes place at the left barrier and the potential 
level in the dot follows that in the right lead. The asym- 
metric behavior of the voltage drop in our model systems 
has analogies with the case of atomic chains between two 
electrodes [33] 
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FIG. 11. Double-barrier potential system B. a) The 
zero-bias voltage effective potential along the symmetry 
a:xis. The energy-zero corresponds to the bottom of en- 
ergy band in an infinite 2D-system with the electron den- 
sity of 0.2 (a*))^. The Fermi level is shown by the dashed 
line, b) The change of 14// due to bias voltage. In the up- 
per panel AVbias ~ 0.03 Ha* (0.36 meV) and lower panel 
AVua, = 0.06 Ha* (0.71 meV). 

The behavior of the potential level in the quantum dot 
is connected to the occupation of the dot resonance state 
and its position relative to the lead Fermi levels. Fig. 12 
shows the local density of states (LDOS) calculated by 
integrating over the quantum dot area. For the zero bias 
voltage, both cases, A and B, have a resonance peak be- 
low the Fermi level. When the bias AVbias is applied the 
potentials and the Fermi levels are shifted by +^ AVbias 
and —^AVbias in the left and right leads, respectively. 
This defines the so-called bias window on the energy axis. 
At small AVbias the value the resonance peak to case B 
moves down in energy. The resonance, which gives a 
large contribution to the charge in the dot is below the 
left Fermi level. The bias induced charge redistribution 
takes place near the left barrier. Thus the potential in 



2000 



5 
S,2000 

§1000 

a- 

g 

at 

5 2000 



Electron energy, <a (meV) 
-0.4 -02 g 02 04 






0.59 0.6 



4000- b) 



2000 



-0.4 



0.61 0.62 0.63 0.64 0.65 0.66 
Electron energy m (Ha ) 



Electron energy, (o (meV) 
-0.2 g 02 0.4 




'4000 



;2000 




0.59 0.6 



0.61 0.62 063 0.64 
Electron energy to (Ha ) 



0.65 0.66 



FIG. 12. LDOS in the region between the barriers shown 
in Fig. 10. a) LDOS for case A with narrow barriers, b) 
LDOS for the case B with wide barriers. The vertical lines 
denote the Fermi level position in the leads. Both in a) and b) 
the uppermost panels correspond to the zero-bias calculation, 
the middle panels to AVuas = 0.03 Ha* whereas the lowest 
panels correspond to AVuas = 0.06 Ha*. 
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The position of the resonance peak relative to the 
Fermi levels has a large effect on the electron transmission 
probability through the double-barrier potential system. 
The current flow is due to the states with energies be- 
tween right and left Fermi-levels i.e. in the bias window. 
When the resonance peak moves into this region there 
is a steep increase in the current. Thereafter the cur- 
rent stays approximately constant as a function of the 
bias voltage. This characteristic behavior of the double 
barrier potential is visible in Fig. 13. Case B with the 
sharper resonance peak has a steeper raise of the current 
than case A. Moreover, the raise occurs at a higher bias 
voltage in case B than in case A. 




0.6 
Bias voltage (mV) 

FIG. 13. Current as a function of the bias voltage for the 
double-barrier potential systems shown Fig. 10. a) Case A 
with the barrier width of lag. b) Case B with the barrier 
width of 1.25 Oq. The zero-bias conductivities of case A and 
B are 0.060 Go and 0.014 6*0, respectively. 



2. Asymmetric barriers 

So far both the potential barriers in the system of 
Fig. 10a have been identical. Inspired by the prospect to 
use non-symmetric molecules as rectifiers [34,35] we have 
studied also double-barrier systems with non-identical 
barriers. The zero-bias conductivities of the cases A and 
B (see Fig. 13 and its caption) are 0.060 Go and 0.014 Go- 
These are of the same order in magnitude as conductiv- 
ities calculated for molecules between electrodes [35]. In 
the next example we have reduced the height of the sec- 
ond barrier in case A by a factor of two in order to create 
an asymmetric system. 

The ensuing current- voltage curve is shown in Fig. 14. 
The curve is asymmetric with respect to the direction 
of the applied bias. The double-barrier system shows a 
clear rectification effect resembling that for asymmetric 
molecular wires [35]. The reason for the rectification ef- 
fect is seen in the LDOS in the quantum dot given in 
Fig. 15. When the bias over the system is zero a res- 
onance peak is below the Fermi level as it was in the 



previous cases A and B. For positive bias voltages (the 
potential is higher in the lower-barrier side) the resonance 
peak moves up in energy and the resonance is emptying 
of electrons. This causes the increase in the conductiv- 
ity. In the case of negative bias voltages (the potential 
is higher in the higher-barrier side) the resonance peak 
follows the Fermi energy of the lower-potential lead. The 
situation is similar to that of system B above at low bias. 
The resonance does not enter the bias window as fast as 
in the case of the positive voltage and the current in- 
creases slowly. 
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FIG. 14. Current-voltage curve for a double-barrier poten- 
tial system with asymmetric barriers. 
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FIG. 15. LDOS for the double-barrier potential system 

with asymmetric potential barriers. The LDOS corresponds 
to the quantum dot region between the barriers. 



V. CONCLUSIONS 

We have developed a computational scheme to model 
two-dimensional nanostructures connected to two semi- 
infinite leads. The electron density and the current 
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are calculated sclf-consistently using the non-equilibrium 
Green's function approch. The single-particle electron 
states are handled within the density-functional theory. 

We have formulated the problem using the finite- 
element approximation. In this approximation the 
boundary conditions are easy to derive and imple- 
ment. We have shown the derivation of the Dirichlet-to- 
Neumann boundary conditions and the discretized forms 
of physical quantities such as the tunneling probability. 

Tests with model potential systems show the numer- 
ical accuracy and its dependence on the finite-element 
mesh chosen. Especially, we show that for efficient ac- 
curate calculation is important to refine the mesh near 
the boundaries between central region and the bound- 
aries. Self-consistent calculations for resonant tunneling 
structures demonstrate the efficiency of the scheme. 

We have treated systems with upto 10 000 degrees 
of freedom. Three-dimensional atomistic systems de- 
scribed by of pseudopotentials would need roughly one 
order of magnitude more degrees of freedom withch is wi- 
hin present-day computational capabilities. The present 
two-dimensional work is an important step in the devel- 
opment towards three-dimensional atomistic modeling of 
non-equilibrium transport in nanoscale devices. 
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